%% boundary_nodes (for Dirichlet)
function bdn = boundary_nodes(mesh, index)
if strcmp(mesh.basis_type, "P2")
    gap_x = 2*mesh.Nx;
    gap_y = 2*mesh.Ny;
else
    gap_x = mesh.Nx;
    gap_y = mesh.Ny;
end
switch index
    case 1 % bottom
        start = 1;
        stop = gap_x + 1;
    case 2 % right
        start = gap_x + 1;
        stop = gap_x + gap_y + 1;
    case 3 % top
        start = gap_x + gap_y + 1;
        stop = 2*gap_x + gap_y + 1;
    case 4 % left
        start = 2*gap_x + gap_y + 1;
        stop = size(mesh.bdn,1);
end
bdn = mesh.bdn(start:stop);
if index == 4
    bdn = [bdn; 1];
end
end